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Abstract 

An upwind solver is presented for the unstructured grid National Combustion Code (NCC). The 
compressible Navier-Stokes equations with time-derivative preconditioning and preconditioned flux- 
difference splitting of the inviscid terms are used. First order derivatives are computed on cell faces and 
used to evaluate the shear stresses and heat fluxes. A new flux limiter uses these same first order 
derivatives in the evaluation of left and right states used in the flux-difference splitting. The k-epsilon 
turbulence equations are solved with the same second-order method. The new solver has been installed in 
a recent version of NCC and the resulting code has been tested successfully in 2D on two laminar cases 
with known solutions and one turbulent case with experimental data. 

Introduction 

The National Combustion Code (NCC) is a state of the art Computational Fluiid Dynamics (CFD) 
program specifically designed for combustion processes 1 . The code employs an unstructured grid, 
massively parallel computing 2 3 , a dynamic wall function with the effect of adverse pressure gradient 4 , a 
low Reynolds number wall treatment 5 , a cubic non-linear k-epsilon turbulence model 6 , a lagrangian liquid 
phase spray model 7 , and stiff laminar chemistry integration. Recently, viscous low-speed 
preconditioning 8 has been added to improve the low-speed convergence of the NCC in viscous regions 
and the ability to handle multiple sets of periodic boundary conditions has also been added. The 
combination of these features is usually not available in other CFD codes and gives the NCC an 
advantage when computing recirculating, turbulent, reacting spray flows. Previously, the NCC has 
undergone extensive validation studies for simple flows 9 , complex flows 10 , NO x emissions prediction 11 
and traditional gas turbine combustor injectors 12 . 

The current NCC solver uses a finite volume discretization with unknowns stored at cell centers. The 
inviscid fluxes are evaluated with average values on the cell faces and a Jameson-type dissipation 
operator is added to avoid odd-even decoupling. For the viscous fluxes, first order spatial derivatives at 
cell centers are obtained by summing cell-face average values and cell-face derivatives are then found by 
averaging adjacent cell -center values. The equations are advanced in pseudo-time by a four stage Runge- 
Kutta scheme in which residual smoothing is applied at each stage. For time-accurate simulations, a dual 
time-stepping approach is employed. 

There are a number of flow regimes in which an upwind solver should produce a more accurate 
solution than the current NCC solver. The most obvious involve transonic flows with shocks. A less 
obvious case involves turbulent shear layers where the decrease in artificial dissipation inherent in the 
upwind method should produce more accurate solutions of the k-epsilon equations. Finally, for reacting 
flows, the decrease in artificial dissipation should again give more accurate solutions for the species 
transport equations. In addition, for these flows as well as others, the use of face-centered derivatives in 
the viscous fluxes should lead to more accurate solutions for all equations. 

In the present work, we introduce an upwind solver for use on the unstructured grids of NCC. 
Preconditioned flux-difference splitting is employed for the inviscid fluxes. First order derivatives are 
computed on cell faces and used to evaluate shear stresses and heat fluxes. A new flux limiter, based on a 
nonlinear average of these same face-centered derivatives, is used in the evaluation of the left and right 
states for flux-difference splitting. Since the test cases in this paper are all two-dimensional non-reacting 
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flows, for simplicity the description of the method is similarly limited. The extension to three- 
dimensional reacting flows is straightforward. 

Governing Equations 

The governing equations are the two-dimensional compressible Navier-Stokes equations with time- 
derivative preconditioning plus the two-equation k-e turbulence model. Let x and y be space dimensions, 
t, p, u, v, h, T, and p are time, pressure, x and y velocities, enthalpy, temperature and density, respectively 
and all variables are dimensional. In addition, for low Mach number flows, it is convenient to subtract off 
the large background pressure p (set by user) by setting 

p = p + p'. (1) 


The equations are written in vector conservation-law form 

P- 1 — + —(E-E )+-(F-F)=H, 
ct ox oy 

Where 

U = {p, pu , pv, e\ pk, pef , 

E = {pu , pu 2 + //, puv, puh t , puk , pusj , 

F = (/tv, puv , pv 2 + p', pvh, , pvk, pvsj , 

E v — (0, ^ xx 1 T xy ’ U ^ XX V t XV — Qx ’ t kx ’ ^ ex ) ’ 

F v — (0, T yx ’ T yy > U T yx V ^ yy ~ Q y ’ T ky ’ ^ ey ) ’ 

H = { 0,0, 0,0, S k ,Sj , 


(2) 


(3) 


P c is the preconditioning matrix, e = ph t — p' is the total energy, and h t = h + \ (u 1 + v 2 ) is the total 
enthalpy. Note that subtraction of the constant background pressure p does not change the basic form of 
the conservation equations. The density p is given by the equation of state for an ideal gas 

p = {p + p')/RT, (4) 

where R is the gas constant for the fluid. The viscous stresses, k and £ viscous terms and the heat flux 
vector are given by 
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here K eff = K + Pr ( 1 p t C p is the effective thermal conductivity, //, p t , K and C p are the fluid 
viscosity, eddy viscosity, thermal conductivity and specific heat, respectively. Also cr k l , (7 1 are model 
constants and Pr f is the turbulent Prandtl number. The turbulence source terms S k , S £ are given by 


S k =P k ~P£, 

where the production term P k is 
\du ; 


s,=cj-p t -c, 2 ^- 


P k =\-pil i U j 


dx 


and C £l , C r2 are additional model constants. The Reynolds stresses and the eddy viscosity are given by 
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Note that we are using only the linear terms in the Reynolds stress model 6 . However, we have found it 
necessary to use the variable C /; from that reference in order to capture the near wall behavior. This is 
given by 
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For turbulent flows the generalized wall function 4 is applied at all solid walls. 

The preconditioning matrix P c is adapted from the work of Turkel 13 . First consider the inviscid 

equations written in terms of the primitive variables V = (p',u,v,S',k,e ) T , where dS ’ = dp’ — c 2 dp, 
c 2 = vRT is the speed of sound and y is the ratio of specific heats. These are written as 
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For this system the simplest form of Turkel’s preconditioner is a diagonal matrix written as 
Pp 1 = Diag(m 2 , 1, 1, 1, 1, l), 


(13) 


where m 2 = j3 2 / c 2 and [J 1 is chosen to bring the system eigenvalues closer together at low Mach 
numbers. The eigenvalues of P p A p are given by 


A. — U — Uq — c 0 , u , u, u, u, u — u 0 + c 0 


(14) 


with 


Cq=uI+J3 2 , u 0 . (15) 

The conservation-law form is then given by P c = MP j: M 1 with M = dU/SV . Finally, following Choi 
and Merkle 8 , we express Eq.(2) in terms of altered primitive variables Q = ( p\ u, v, h, k, sf 

(16) 


Pq — + —{e-e v ) + — (f-f v )=h 

Q dt dx y vJ dy y vJ 


where P Q 1 = MP p 1 N and N = GV / GQ . 

For most flows it is sufficient to set j3 2 equal to the square of the local flow speed together with 
lower and upper limits. However, for viscous flows near walls, convergence is improved by introducing 
the viscous limit of Choi and Merkle 8 . Hence we take 

P 2 = min [max [fiy , u 2 +v 2 , p 2 ), c 2 ]. 


( 17 ) 


where p 2 is a global user specified lower limit and 
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( 18 ) 


2 _ a(a- l)v,; 

~ \ + {a-\)vl / c 2 

Here a = /?CV /Re cd; , RCV is the ratio of CFL to von Neumann numbers (user specified), 

R e cd , = pv n 8 n j j.L e(( is the local cell Reynolds number, v n is the velocity normal to the cell face and 8 n 
is the spacing between adjacent cell centers. 

Discrete Formulation 


The governing equations are solved on an unstructured grid (triangles or quads in 2D) with 
unknowns stored at cell centers. Integration of Eq. (16) over a cell together with explicit time differencing 
gives 

P Q l (AA/ At )AQ = -R" , (19) 

where AA is the cell area, At is the time step, A Q = (9" —Q", and the area- weighted residual R is 
given by 

R , -(r-r.lAj-tfM • < 2 «> 

where the summation is over all faces j of cell i , while dy tj , dx jj are the changes in y,x respectively, 
along face j. 

The cell-face inviscid fluxes E and F in Eq. (20) are approximated by preconditioned flux- 
difference splitting 14 

e = r[e{u l )+e{u r ^-\p; 1 \p c a c \u r -u l ), ( 21 ) 

where U L and U R are, respectively, left and right states, A r = dE/dU , A c denotes evaluation at an 
average state, and P r 1 1 P c A c | = MP ^P p A p | M 1 . Here \P p A p | is computed using the eigenvectors and 
absolute eigenvalues of P p A p . Eq.(21) is rewritten in terms of the altered primitive variables Q 

E = p[e(q l )+e(q r \-\A°(q r -Q l ), (22) 

where A“ = MP~ l \P p A p | N. The left and right states Q L , Q R are obtained by extrapolating from the 

adjacent cell -centers. Note this linear extrapolation leads to a second order approximation for the fluxes. 

Let i\, i2 denote the centers of two adjacent cells and fc denote the centroid of the cell face separating 
these two cells. Then 

e‘=e„ + f^ 




In the present implementation the average state Q is taken as a simple average of Q L and Q R . The flux 
F is approximated by expressions similar to Eqs. (21) and (22) with A c replaced by B = dF/dU and 
A a replaced by B" = MP p l \P p B p \N . 

The cell-center derivatives in Eq. (23) are obtained by a weighted average process that is a 
modification of that used by Huynh 15 for triangular cells. Let ij represent face j of cell i and set 
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where these expressions are applied to each of the variables in Q . Note the resulting cell-center 
derivative is biased toward the least steep of the surrounding cell-face slopes and thus serves as a flux 
limiter. Also, this formulation works for any cell shape, in two or three dimensions. 

The cell-face derivatives are obtained by integrating over a volume surrounding the face. Let i\, i 2 
denote the centers of two adjacent cells and kl, k 2 the ends of the face separating these two cells. 
Averaging the derivatives over this volume gives 

it] =i <M(e,2 -QnW -(e« -Q u Wi 

\ OX J f 


[(e„ - e„ ]dx“ - fo, 2 - e„ W. ■ 


where 


dx‘ = x i2 - x n , dy' = y i2 - y a , dx k = x k2 - x kl , dy k = y k2 - y kl , 
dA f = 4 ( dx'dy k - dy'dx k ) . 

Since Q k 1 and Q k2 are not stored on the grid, their difference has to be approximated. First we obtain 
simple approximations to the cell-center derivatives by integrating over the cell volumes, 


= cia : 1 Y,Q/iy J . ^ ] = -<«-■ 'Lay * 1 , 

8x). “ \Sy), j 


where the sum is over all the cell faces and the face values Q j are taken as the average of the adjacent 
cell-center values. Then we approximate Q k2 —Q kl as 


. m ( dQ) \ k , ae [ae . » 




Eqs. (26) together with Eq. (29) are used to evaluate the cell-face derivatives. These are used not only for 
the weighted averages of Eq. (22) but also for the cell-face shear stresses and heat fluxes of Eq. (5). 

Time Marching Solution 

Eq.(19) for Q is advanced in pseudo-time by a four stage Runge-Kutta scheme. The viscous stresses 
and heat transfer in Eq. (5) together with the numerical dissipation in Eq. (22) are updated in the first two 
stages only. For the turbulence variables k and £ , the updated values are under-relaxed at each Runge- 
Kutta stage. Finally, residual smoothing is applied at each stage using two sweeps of an under-relaxed 
Jacobi iteration. 
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For the steady flow problems considered here, convergence is based on the root-mean-square of the 
X and y momentum residuals, res tot . The solution is considered converged when res tot / res) ot < goal 

for 100 consecutive iterations, where res) ot is the value after 1 iteration and goal = 5 x 10 4 for the two 

laminar cases and goal = 10 5 for the turbulent case. For all cases the CFL number is set to 4 while the 
Von Neumann number is set to 1 for the laminar cases and 0.5 for the turbulent case. 

Results and Discussion 

Three test cases have been chosen to demonstrate the accuracy of the upwind method. For the two 
laminar cases, existing numerical results are used to check the current solutions. For the turbulent case the 
present results are compared with experimental data. Note, for the laminar cases, we assume air at 1 atm. 
and 300 K with a reference length D of 1 m. The reference velocity u ref is then chosen to give the 

desired Reynolds number. Re = pu ref D j ju . 

Lid-driven cavity 

The lid-driven cavity is a common benchmark for recirculating flows. The second-order 
streamfunction-vorticity results of Ghia et al 16 ., computed on a uniform 256 by 256 grid, are generally 
accepted as standard. Isothermal flow is set up in a square cavity, side 1 , with a top that moves to the right 
at constant speed u ref . The current work uses a triangular grid with 128 cells along each wall and grid 

spacing along the wall varying from 0.0005 in the corners to 0.02 in the middle. The resulting grid, with 
10,514 triangular cells, is shown in Fig. 1. 

The case Re = 1000 is obtained by setting u ref = 0.015936 m/sec. A contour plot for uju ref , shown 

in Fig. 2, gives the general structure of the flow. Note the two recirculating regions in the bottom corners 
of the flow. Convergence to the desired level is shown in Fig. 3. Note that this took 6683 iterations. 
Profiles of u/u ref on the vertical centerline and vju^ on the horizontal centerline are compared with 

the Ghia et al. results in Figs. 4 and 5, respectively. Considering the sparseness of the grid, the 
comparison is excellent. 

Laminar backward-facing step 

We consider isothermal laminar flow over a backward-facing step with a step height of half the 
downstream channel height D and a parabolic inlet velocity u(y) = 24y(0.5 — y)u ref f° r 0 < y < 0.5 

specified at the step. The outflow boundary is taken at 15 channel heights downstream of the step. 
Gartling 17 set up this problem as a test for outflow boundary conditions. Numerical results were obtained 
for Re = 800 with the outflow boundary at 30 channel heights and used a Galerkin-based finite element 
method with 800 x 40 biquadratic velocity elements and linear discontinuous pressure elements. The 
current work uses a 240 by 64 quad grid with stretching in y only. The grid spacing is set at 0.001 at the 
top and bottom walls and 0.03 at the centerline. 

The case Re = 800 is given by setting u ref = 0.012749 m/sec. A contour plot for uju ref , shown 

in Fig. 6, gives the structure of the flow. Note the separation bubble on the upper wall at x « 7 . 
Convergence to the desired level is shown in Fig. 7. Profiles of ulu ref at x — 7 and x — 15 are 

compared with Gartling’s results in Figs. 8 and 9, respectively. Again the results are quite good. Note that 
X — 15 is the exit boundary for the current computation. 
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Turbulent backward-facing step 


Kim, Kline and Johnston 18 ' 19 (KKJ) studied the turbulent flow over a backward-facing step with two 
main objectives. The first was to characterize the flow in the separated shear layer, the reattachment 
region, and the developing boundary layer downstream of reattachment. The second was to obtain 
experimental data that would be useful in the development of computational models. In KKJ’s experiment 
the step height h was 3.81 cm, the channel width upstream of the step was 2 h , and boundary layer trips 
were placed on each side wall at a distance of 8 h upstream of the step. On the centerline of the channel, a 
distance of 4 h upstream of the step (termed the reference point), the reference speed u ref was 18.2 m/s 
and the dynamic pressure was 20.3 mm of water. 

In the current work the computational domain is taken from 8 h upstream to 30 h downstream of the 
step. We use a stretched quad grid with 30 by 32 cells (x, y) upstream and 50 by 64 cells downstream. 
The y grid spacing is set at 0.002 h along both side walls. The x grid spacing varies from 0.1/7 at the 
upstream boundary (— 8 h) to OAh at the step, remains constant to x — lOh , and then increases gradually 
to 1.312/t at the exit (3O/7). The resulting grid, from x/h = -4 to 20, is shown in Fig. 10. The 
computation is started at — 8/7 to allow the turbulent boundary layers to grow to suitable thickness before 
reaching the reference point. At the upstream boundary we assume a uniform inflow with a speed u m , 
turbulence level of 5%, and a mixing length of 0.001m. The inlet speed was varied until the solution at the 
reference point agreed with the experimental value u ref — 18.2 m/sec. This required a value of u in = 

17.5 m/sec. Note this value was not very sensitive to the inlet values of turbulence level and mixing 
length. A contour plot of u / u ref , shown in Fig. 1 1 , gives the general structure of the flow and 
convergence to the desired level is shown in Fig. 12. The computed reattachment length, as determined 
from the change in sign of the wall shear stress, was x R /h= 6.75. This is well within the experimental 
range of 7 + 1 reported by KKJ. 

The computed wall-static pressure on the step side of the channel (plotted as the pressure coefficient 
Cp ) is compared with experimental values in Fig. 13. Here Cp is given by 


_ P-Pref 

\Pulef ' 


(30) 


We determine p ref by forcing the computed pressure to go through the first data point downstream of the 
step. This comparison is quite good although the computed pressure rise occurs a little early. Profiles of 
u/u re f at four locations downstream of the step are compared with the experimental values in Figs. 14 to 

17. Again the comparison is quite good except in the reverse flow portion of the separation zone. Here 
KKJ note that the flow was unsteady and the measurements may not be accurate. 


Conclusions 


An upwind solution method for the unstructured grid National Combustion Code has been presented. 
The solver uses preconditioned flux-difference splitting for the inviscid terms and face-centered first 
derivatives for the shear stresses and heat fluxes. A new flux limiter uses these same first order 
derivatives in the evaluation of left and right states used in the flux-difference splitting. The same 
methodology is used for the k-epsilon turbulence equations. The solver has been tested in 2D on two 
laminar and one turbulent flow and shown robust convergence in all cases. The laminar results appear 
quite accurate when compared with known computational solutions. The turbulent results are in good 
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agreement with the experimental data. The second order solution of the k-epsilon equations shows 
promise for the future analysis of turbulent reacting flows. 
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Fig. 1. Grid for lid-driven cavity with 10,514 triangular cells. 
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Fig. 2. Contours of u/u re f for the lid-driven cavity at Re = 1000. 

10 






Fig. 5. Comparison of v/u re f on the horizontal centerline with the results of Ghia et al for Re = 1000. 
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Fig. 6. Contours of u/u re f for the laminar backward-facing step at Re = 800, 
(where the grid portion from x = 10 to 15 is not shown). 


12 


10 c 


10' 1 v 


w 

CD 


° ? 

(/> 10 2 


CD 


1 0 3 1 


0 


Fig. 7. Conve 
0.5 
0.4 
0.3 
0.2 
0.1 

> 0 

- 0.1 

- 0.2 

-0.3 

-0.4 

-0.5 


Fig. 8. 





Fig. 9. Comparison of u/u re f at x — 15 with Gartling’s results at Re = 800. 
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Fig. 10. Grid for turbulent backward -facing step 
(full grid goes from x/h = -8 to 30) 
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Fig. 1 1 . Contours of u/u re f for turbulent backward-facing step 



Fig. 12. Convergence of resto/rest,,/ for the turbulent backward-facing step. 
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